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We present a detailed theoretical study of the transfer of electronic excitation en- 
ergy through the Fenna-Matthews-Olson (FMO) pigment-protein complex, using the 
new developed modified scaled hierarchical approach [Shi Q. et al, J Chem Phys 2009, 
130, 084105]. We show that this approach is computationally more efficient than the 
original hierarchical approach. The modified approach reduces the truncation levels 
of the auxiliary density operators and the correlation function. We provide a sys- 
tematic study of how the number of auxiliary density operators and the higher-order 
correlation functions affect the exciton dynamics. The time scales of the coherent 
beating are consistent with experimental observations. Furthermore, our theoretical 
results exhibit population beating at physiological temperature. Additionally, the 
method does not require a low-temperature correction to obtain the correct thermal 
equilibrium at long times. 



I. INTRODUCTION 

In the initial step of photosynthesis, light is captured by protein-bound pigments that 
are part of light-harvesting antenna complexes. The excitation energy is transferred with a 
near-unity quantum yield to reaction centers, where it is converted to chemical energy. The 
underlying molecular mechanisms responsible for the near-unity quantum yield are far from 
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being understood. In this regard, an extensively studied and relevant system is the Fenna- 
Matthews-Olson (FMO) pigment-protein complex of green-sulphur-bacteria, which acts as a 
mediator of excitation energy between the outer antenna system, i.e., the chlorosomes, and 
the reaction center lfl. 

Experimentally, Savikhin et al. [2] observed quantum beating in the FMO complex 
using the fluorescence anisotropy technique. More recently, Engel et al. [3] employed two- 
dimensional electronic spectroscopy to observe long-lasting quantum beats that provides 
direct evidence for survival of long-lived electronic coherence for hundreds of femtoseconds. 
Aspuru-Guzik et al. investigated the effects of quantum coherence and the fluctuat- 

ing environment, using the Lindblad formalism, on the enhancement of the photosynthetic 
electronic energy transfer efficiency from the perspective of a quantum walk. In Rebentrost 
et al. jsL a method to quantify the role of quantum coherence was introduced. Ishizaki 
et al. [7H9| employed the hierarchical equation of motion (HEOM) approach expansion in 
order to address the robustness and the role of the quantum coherence under physiological 
conditions. Their results reveal that quantum wave-like motion persists for several hundred 
femtoseconds even at physiological temperature T = 300K. Very recently, a very large-scale 
calculation of energy transfer between chromophore rings of purple bacteria was carried out 



using the HEOM approach [10fl. Meanwhile, Whaley et al. 11] and Caruso et al. 12 1 
discussed quantum entanglement in photosynthetic harvesting complexes and clarified the 
connection between coherence and entanglement. They showed that the FMO complex ex- 
hibits bipartite entanglement between dimerized chromophores. The subject continues to be 
of great interest with a large number of publications discussing electronic energy transfer in 
photosynthetic complexes, in particular, the issue of quantum speed up in the FMO complex 
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19J. 



In this paper, we present a detailed theoretical study of the transfer of excitation energy 
towards the reaction center through the Fenna-Matthews-Olson (FMO) pigment-protein 
complex using a modified scaled hierarchy equation of motion approach, which was devel- 
oped by Shi and coworkers recently [20 ] ■ This approach guarantees that the auxiliary density 
operators decay to zero at high truncation level. Furthermore, it provides a considerable 
computational speedup over the original hierarchical approach We will show that the 
scaled hierarchical approach can reduce the truncation level of both auxiliary density op- 
erators and the correlation functions compared to the classical approach. The time scales 
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of the coherent beating are consistent with experimental observations. Furthermore, our 
results show that the population beating persists at physiological temperature. 

II. THEORY 

Excitonic energy transfer within photosynthetic proteins -such as the FMO complex- op- 
erates in a demanding parameter regime where a small perturbative quantity is not available. 
It is thus a challenge to find accurate and efficient methods for the simulation of the quan- 
tum dynamics. A number of approximate methods have been developed |16|: they include 
the semi-classical Forster theory, standard Redfield theory, modified Redfield theories, the 



modified Lindblad formalism and hierarchical equations of motion, among others |4j, |2l|, |22|. 
In this study, we utilize the recent hierarchical Liouville space propagator method de- 

energy transfer in the FMO complex. 
The original method is based on a reformulation of the original hierarchical quantum mas- 
ter equation and the incorporation of a filtering algorithm that automatically truncates the 
hierarchy with a preselected tolerance. They showed how this method significantly reduces 
the number of auxiliary density operators used to calculate electron transfer dynamics in a 
spin-boson model and the absorption spectra of an excitonic dimer 



20]. 



The structure of FMO complex was first analyzed by Fenna and Matthews in 1975 23]. 
It consists of a trimer, formed by three identical monomers. Each monomer contains seven 
bacteriochlorophyll a (BChl a) molecules (or seven "sites".) Biologically, the FMO complex 
acts as a molecular energy wire, transferring excitation energy from the chlorosome structure, 
where light is captured, to the reaction center (RC). There is substantial evidence that the 
FMO complex is oriented such that sites 1 and 6 are close to the baseplate protein and 
sites 3 and 4 are close to the RC complex and thus define the target region for the exciton 
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23H25]. The detailed structure is shown in Fig. [TJ It should be mentioned that the 
presence of the eighth BChl a molecules per monomer has been proposed by Ben-Shem et 
al in 2004 26]. This has been verified experimentally recently by Tronrud and coworkers 
27] . It has been suggested that the eighth BChl a molecule acts as a gateway site from the 
reaction center to the 27 chromophores in the trimer [27]. 
The total Hamiltonian of the quantum system is given by: 
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n — Us + Tib + n.sBi 



(1) 



where Tis, "Hb, and T-Lsb are the Hamiltonian of the system, the environment, and the 
system-environment coupling, respectively. Here, we consider each site as a two-level sys- 
tem of ground state and excited state. The system Hamiltonian, Hs, which describes the 
electronic states of the pigments can be expressed as: 



where |j) denotes the state with only the j-th site is in its excited state and all other sites 
are in their ground state. Sj represents the site energy of the j-th site which is defined 
as the optical transition energy at the equilibrium configuration of environmental phonons 
associated with the ground state. N is the number of pigments or sites. Jj^ is the electronic 
coupling between site j and k. The parameters for this Hamiltonian are taken from the 
paper of Adolphs and Renger [28]. In their work, two independent methods were used to 
obtain the site energies of the seven BChl a molecules of the monomeric subunits of the FM O 
complex. In the first method, the site energies are used as parameters that were optimized 
by a genetic algorithm in the fit of the optical spectra. In the second one, the site energies 
are obtained directly by electrochromic shift calculations j^sj. 




For the Hamiltonian of the environment, a harmonic oscillator model is applied. 

Furthermore, it is assumed that the electronic excitation on each site couples to its own 
bath independently: 



where Njb is the number of different harmonic modes coupled to the j-th site. Here, m,j£, 
Pj^ and Xj^ are mass, frequency, momentum, and position operator of the harmonic bath 
modes. The coupling term, Hsb, which is responsible for fluctuations in the site energies by 
the phonon dynamics, can be expressed as: 
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(2) 
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where the bath part is Fj = Cj^-Xj^ and the represent the system-environment coupling 
constant for the j-th site and £-th phonon mode. The projection operator Vj = 
describes the system part of the interaction. 

At time t = 0, we assume that the system and the environment are decoupled, i.e. 
Ptot (0) = p (0) <S> Pb (0). Additionally, the environment is in the Boltzmann equilibrium 
state, Pb (0) = e ~ pHB /Tr B [<e~P HB ], where = 1 /k B T. The time evolution of the system density 
matrix, p(t), can be calculated by tracing out the environment degrees of freedom: 



p (t) = Tr B [p tot (*)] = Tr B p tot (0) <e^] . ( 5 ) 



The bath is described by its correlation functions, Cj (t), which are defined as 1C| |29|, |30|, 
Cj (t) = Trs[Fj (t) Fj (0) ps], where the Langevin force, Fj (t), is given in the interaction 
picture: Fj(t) = <& lHBt ' h Fj <e~^ s *' R . For the phonon bath, the correlation function can be 
written as Cj (t) = ^ du ■ Jj (u) ■ 1 ^ , where Jj (u) is the spectral density for the 
j-th site: 

jj h = £ o^r 5 ^-^- (6) 

$ -'".it ■ U K 

To proceed, we use the Drude spectral density, which corresponds to an overdamped Brow- 
nian oscillator model. Furthermore, we assume that the system-environment coupling is the 
same for all sites, Jj (u) = J (u), V js. The Drude spectral density is defined as: 

J H=^^-r (7) 

We introduced r] = ^ -which is dependent on the reorganization energy A- and the Drude 
decay constant, 7. Under this spectral density, the correlation function Cj (t) takes the form: 

00 

C j (t>Q)=^2c k -e-^ i , (8) 



fc=0 



with the Matsuraba frequencies vq = 7 and V/. = ^ for k ^ 1. The constants are given 
by: 
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/or k ^ 1. 



Now, we are in the position to write down the HEOM for the reduced density operator jlO|. 
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dt' 
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(9) 
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where n denotes the set of nonnegative integers n 



{ni,n 2 , • • • ,n7v} 



{{"io, n X i 



n iK}, • ■ • > {^Mb ^tvi, ' ' ' , n Af_ft:}}- n ^ refers to the change of the number 



rijk to rijk ± 1 in the global index n. The sum of rijk is called tier (M c ), M c = Yljk n J k - ^ n 
particular, p — P{{o,o,-},-,{o,o,-}} is the system's reduced density operator (RDO) and all 
others are auxiliary density operators (ADOs). Although the RDO is the most important 
operator, the ADOs contain corrections to the system-environment interaction; these arised 
from the non-equilibrium treatment of the bath. 

Here, we assume that both p Q and Vj have the order of one. When the tier of p n at 



tier M c (M c = ££ k n jk )> ^e amplitude o: 
following the standard approach (Eq. [9]) 



P 
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is proportional to 



■-K 



32], which indicates the amplitude of p n is 



related to both Ck and n^. The rijk is decided by the truncation level, while the c& is related 
to the correlation function. The correlation function is derived from the system-environment 
correlation. In other words, the amplitude of p n is dependent on the system-environment 
coupling. Under the intermediate-to-strong system-environment coupling, the amplitude of 
p n can not be guaranteed to be small even at high truncation level. It goes to the opposite 
direction as we expected as we always expect more accurate results at high truncation level. 
Fortunately, Shi and coworkers developed a new approach in which one is able to rescale the 
original ADOs which can be used for overcoming this issue [20]. They scaled the original 
operator as: 
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Pn (t) 



Pn (*) , 



(10) 



After the scaling, the \p n \ has the order of \\ k . ^/^^Jn^i . It can make sure that |p n | 
decays to zero at higher hierarchical truncation level. 

Since the number of contributing terms to the correlation function, Eq. [HJ and ADOs 
are infinite, the computation of Eq. is -in general- impossible. In order to overcome this 
problem, a truncation scheme for both the correlation function and ADOs is applied. We set 
the truncation level for the correlation function (Matsuraba frequency and constant c^) at 
level K, while the cutoff for the tier of ADOs is M c - With the Ishizaki-Tanimura truncating 



scheme 
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32], Eq. [9] for the scaled density operator becomes: 



d 



N K 

- ( iC s + ^ ^ n jk v k J p, 



A? 



" z V + !) k*| [V i5 p n + 

Af oo 

-E E ^ • [v„ [v„ pj 

-i E E (cfcV,- p n - - 4p n - V,-) . 



3=1 fc=0 

We use Eq. [TT] to simulate the exciton dynamics of the FMO complex. 



(11) 



III. RESULTS AND DISCUSSION 



For the numerical analysis, we used the same Hamiltonian as in Refs. |7J, [28|, the same 
reorganization energy, Aj = A = 35 cm -1 , and Drude decay constant, 7 J ~ 1 = 7 -1 = 50 fs, 
of Ref. j3, 33|, As we mentioned before, sites 1 and 6 are both connected to the LHC It is 
possible that sites 1, 6 or both are excited. For this reason, three different initial conditions 
are employed, 1 1) (site 1 is excited), |6) (site 6 is excited) and (|1) + |6)) (the superposition 
of excited sites 1 and 6). 

Calibration— We compare the scaled approach to the original HEOM approach and inves- 
tigate the critical choice of the truncation levels of ADOs (AQ and the correlation functions 
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(K). The original HEOM approach is given in Eq. ([9]). In Fig. [21 we depict the population 
of sites 1, 2 and 3 for different N c in the scaled approach and for N c = 4 in the original 
approach at temperatures T = 77 K and T = 300 K, setting K = 0. One obtains a large 
difference between the two approaches at T = 77 K, see the dotted lines in the figure. Un- 
der the original HEOM method, the population of site 2 goes below after 750 fs, which is 
unphysical. However, the population of each site under the scaled HEOM approach behaves 
reasonably even at M c = 1. This shows that the scaled HEOM approach can result in better 
simulation results at less computational costs. The difference between the two approaches 
originates from the truncation level of the correlation function, which is due to the coupling 
between the system and the environment. 

For the scaled approach itself, the difference between different truncation levels (A/" c ) is 
modest at both temperatures. It can be seen that there is only minimal difference among 
three Af c values at T = 77 K. Beyond 1500 fs, the difference of the population evolution 
for site 3 becomes larger between the case of Af c = 1 and M c = 2 and 4. The population 
evolution of all the sites is exactly the same for M c = 2 and 4. As a result, M c > 2 is the 
sufficient truncation level for the ADOs at T = 77 K. At room temperature, T = 300 K, 
the variation at the dynamics of the populations as a function of N c is more apparent. The 
population evolution of all sites for Af c = 1 is not as smooth as in the other two situations. 
For the cases of M c = 2 and 4, there is a slight difference in the population beatings which 
occur between 200 fs and 300 fs. Although this is not a substantial difference, we believe 
M c = 4 is good compromise between efficiency and accuracy. 

For the truncation level of the correlation function (K) , the simulation for the seven sites 
is computationally unwieldy. Therefore, we truncated the system to test the correlation 
truncation level using a three-site model (sites 1, 2 and 3) with three different values of 
K, being K = 0, 1 and 2. The results show that truncation level K = is enough for 
both T = 77 K and 300 K. A similar result was also found in [20], where non-zero K was 
shown to be significant in the dynamics only at rather long time scales. In the following 
computations, we choose M c = 4 for both temperatures as our reference. At this point, 
we would like to emphasize the numerical efficiency of the scaled HEOM approach. The 
original HEOM approach requires a truncation level as high as M c = 12 to get converged 
results jj]. However, we only require J\f c = 2 for T = 77 K and J\f c = 4 for T = 300 K, which 
is a significant resource reduction. On a standard desktop computer, a simulation of the 
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time-evolution for 2.5 ps takes about 7 minutes for the case of M c = 2 and about 1.5 hours 
for the case of M c = 4. 

Coherent beatings at cryogenic and room temperatures— Now, we investigate the cryogenic 
temperature T = 77 K in more detail. This being the temperature of the first experiment by 
Engel et al. [3] which shows coherent phase evolution of the FMO complex from time t = 
to roughly t = 660 fs. The results are presented in Fig. [3l On the left panel, we show the 
results of simulation for the system Hamiltonian only. The the right panel one observes that 
the quantum beating between certain sites clearly persists in the short time dynamics of the 
full FMO complex. For the simulated initial conditions, the population beatings can last for 
lundreds of femtoseconds; this time scale is in agreement with the experimental observation 
sjl . The population beating for all three different initial conditions can last around 650 fs. 
In Fig. |3]^a), the initial state is localized at site 1. The system exhibits coherent beatings 
between the strongly coupled sites 1 and 2, accompanied by relatively slow relaxation to 
sites 3 and 4. The change of population of all other sites is weak. In Fig. [3]^b), where 
the initial state is localized at site 6, the population relaxes faster. For t < 400 fs, there is 
population beating between the strongly coupled sites 6 and 5, accompanied by relaxation to 
the intermediate sites 4, 5 and 7. From these sites, the population is fed into the low-energy 
sites 3 and 4. The population of site 6 almost vanishes at t = 800 fs, while for the previous 
initial condition, Fig. [3]^a), the population of site 1 is roughly 0.5 at that time. The exciton 
migration pathways and time scales are in accordance with previous work 28 1. Finally, 

Fig. [3](c) represents the superposition of site 1 and 6. The time evolution of this case is 
the combination of the single site excited cases. That is the population evolution on site 1 
follows the pathway of initially single site 1 excited, while the pathways for the population 
on site 6 is the same as the single site 6 excited case. 

In order to investigate the excitation transfer beyond the initial beating region, we ex- 
tended the simulation to a longer time (~ 2500 fs). The result are shown in Fig. [5j To check 
whether the entire system converges to the thermal equilibrium, we obtained all eigenvalues 
of the system Hamiltonian and calculate the probability of each eigenstates under tempera- 
ture T based on the Boltzmann distribution. Subsequently, we transformed the population 
from the eigenstate representation to site representation and obtained the population of each 
site at thermal equilibrium. At T = 77K, the population of site 3 is 0.69 and that of site 
4 is 0.22. The population of all other sites is smaller than 0.03. For the case of having the 
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initial excitation start in site 6, the thermal equilibrium is reached the end of 2ps, while 
for the case of the other two initial conditions, they are still on their way to the thermal 
equilibrium at the end of 2.5ps. Our simulation shows that the system reaches thermal 
equilibrium at ~ 7 ps for the case in which site 1 was initially excited and the time for the 
initial superposition of site 1 and 6 is around 6 ps. 

Recent experiment studied the excitation dynamics of the FMO complex at room tem- 
perature {34)]. To investigate quantum coherence effects under physiological conditions, we 
simulate the dynamics at the temperature T = 300 K. We choose three different values 
for 7 _1 in our calculation: 50 fs, 100 fs, and 166 fs J7, 9|. Following the same procedure as 
before, we consider three different initial conditions with the reorganization energy A = 35 
cm -1 . We choose the truncation level at Af c = 4. The calculation results are shown in Fig. 
6. 

The main difference between the case of 300 K and that of 77 K is the time scale of 
the persistence of population beating. The coherent beating lasts only 400 fs at room 
temperature whereas it lasts much longer at T = 77 K. It is also found that the smaller 7 
is, the longer the population beating can last. When 7" 1 = 166 fs, the population beating 
time can last almost 700 fs. The main pathways for all cases are the same as the pathways 
at low temperature. Furthermore, the time evolution of the population for each site also 
converges to thermal equilibrium. However, the entire system reaches thermal equilibrium 
considerably sooner at room temperature. For example, the system initialized at site 6 
reaches equilibrium at 1.5 ps when T = 300 K, compared to 2 ps at 77 K. 

Behavior of the auxiliary density operators— In order to further investigate the effects of 
the ADOs and their role in the modified HEOM scheme, we examined the magnitude of their 
population elements and found that the majority of them are close to 0. However, there are 
some non-zero ADO elements during the time evolution. We plot their time dependence in 
Fig. HJ For the case where site 1 is initially excited, the simulation shows that the most 
important ADO elements are (l|pnoooooo|l) with n — 1, 2, 3, and 4. While for the site 6 
initially excited case, the important ADOs are (6|pooooomo|6) with m = 1, 2, 3, and 4. From 
the image (Fig. Hj), it can be found that the amplitude of the ADOs decays rather quickly 
as the level of truncation increases. Conversely, the ADO populations are related to the 
amplitude of the site population in RDO. When the population goes up, the corresponding 
population in ADOs also increases. For odd truncation levels (pioooooo, P3000000, P0000010 and 
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Pooooo3o) ; the ADOs yield negative population. The results are indicative of the fact that 
the scaled HEOM approach reduces the amplitude of ADOs as truncation level increases. 
Interestingly, there are no negative population elements of the density matrix when the 
truncation level of the correlation function is K = at cryogenic temperature. This is 
in contrast to the original hierarchical approach [7|, in which some populations become 
negative. 

Energy transfer pathways— We briefly comment on the exciton transfer pathways. The 
pathways are determined by the system Hamiltonian rather than by the system-environment 
coupling or the environment From Fig. [3l[5l|6] and [7], we can find that the frequency 

of site population oscillation is independent of the temperature and different bath relaxation. 
For example, in the pathway site 1 ^ 2 — > 3&4 the main beatings between site 1 and 2 
are caused by several features. The energy barrier between site 1&2 (Ae 12 = —120 cm -1 ) is 
smaller than that of site 1&6 (Ae 16 = —220 cm" 1 ) and the coupling of sites 1&2 is stronger. 
The population oscillation between site 3&4 is due to the similar site energies and the strong 
coupling between them. In the real biological system, the RC is close to site 3 and 4. When 
the exciton transfers to these sites, it moves to the RC directly and cannot return to the 
system. Other pathways start from site 6, i.e. site 6^7 — > 3&4, site 6^5 — > 3&4 
and site 6 — > 3&4. Although the population distribution is the same for all pathways at 
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For 



long times, the excitation transfer time is shorter for an initial excitation at site 6 
the pathway of site 1, the site energy of site 2 is bigger than that of site 1. It is hard for 
the excitation to move from site 1 to site 2 and that is why the wave-like evolution lasts for 
a longer time. However, in the population pathway that starts from site 6, the excitation 
flows from the higher- to lower- energy sites all the time. This reduces the transfer time and 
the system reaches thermal equilibrium faster. Comparing the two different temperatures 
in terms of the excitation transfer pathways, we note that the influence of the thermal bath 
is much stronger at room temperature than at low temperature. The bath at 300 K helps 
the system to transfer the excitation more efficiently by reducing the quantum beating, thus 
speeding up the overall transfer times to the reaction center. 
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IV. CONCLUSION 



In summary, we have examined the full dynamics of the transfer of excitation energy 
towards the reaction center through the Fenna-Matthews-Olson (FMO) pigment-protein 
complex, employing the modified scaled hierarchical approach recently developed by Shi et 
al. 20|. The scaled HEOM approach not only reduces the cutoff for the tier of auxiliary 
density operators, but also decreases the truncation level of the correlation function, which 
makes it more efficient compared to the original HEOM approach. We have shown that 
a tier cutoff of N c = 4 and a correlation function cutoff of K = optimizes simulation 
efficiency and accuracy for the parameter regime of the FMO complex. Furthermore, our 
theoretical results show that the population beating can last as long as 650 fs under cryogenic 
temperature (77 K). When the temperature is 300 K, the beating time can vary from 400 fs 
to 700 fs, depending on the environment parameters. Our simulation result is in accord with 
the conclusion of Ishizaki et al [7j. The improved computational performance of our scaled 
HEOM approach will be especially useful in theoretical studies of transport measures such 
as: efficiency; transfer time; and other properties, such as entanglement. Moreover, this 
efficient approach also provides us with the potential to couple other effects into our current 
system. Under the current model, only the thermal effect is fully considered; however, there 
exist many other effects in the real biological system such as: dipole-dipole interaction; the 
different phonon environment for each site; and slow structure changes of the FMO complex. 
It will be our future task to build a model with these features and examine the time evolution 
of entanglement and related quantum information measures 35 
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Figure 1: Sketch of the energy flow in the process of photosynthesis. The energy is captured by the 
light-harvesting complexes (LHC) and transferred to the reaction center (RC). The FMO complex 
is the link between LHC and RC and it operates as a "wire" during the energy transfer process. 



Using the convention for numbering the BChl a molecules (sites) of FMO complex as in ref. 23| , 
site 1 and 6 are close to the LHC and site 3 and 4 are close to the RC. In our theoretical description, 
the excited states of the BChl a molecules of the FMO complex are considered as the "system" and 
all the other relevant degrees of freedom are referred to as the "environment". 
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Figure 2: The population evolution of site 1, 2 and 3 under different cutoffs for the tier of auxiliary 
density operator (AQ and different HEOM approaches. The solid lines represent the population 
evolution at three different truncation levels N c = 1, 2 and 4 of the scaled HEOM approach. The 
short dot lines show the time evolution of N c = 4 at the original HEOM approach. Site 1 is 
initially excited and the reorganization energy and Drude decay constant are Xj = X = 35 cnr 1 
and 7" 1 = 7 -1 = 50 fs, respectively. The dynamics are shown at cryogenic temperature T = 77K 



(upper panel) and at physiological temperature T = 300K (lower panel). 
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Figure 3: The population evolution of each site at cryogenic temperature, T = 77 K. The left panel 
shows the dynamics for the system alone and the right includes the effects of the environment. 
The reorganization energy is Xj = A = 35 cm -1 , while the value of Drude decay constant is 
T^ 1 = 7 _1 = 50 fs. The initial conditions are site 1 excited (a), Site 6 excited (b) and the 
superposition of site 1 & 6 (c). 
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Figure 4: The population evolution of RDO and ADOs for the case of initial excitation at sites 1 
and 6 respectively. The first panel shows the time evolution of the ADO (l|pnOOOQOo|l) elements with 
n = 0, 1, 2, 3, and 4. The second panel shows the time evolution of the (6|/9ooooOmo|6) elements 
with levels of truncation m = 0, 1, 2, 3, and 4. 
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Figure 5: Long time-dynamics of the population at each site for T = 77 K, where (a) , (b) and (c) 
are corresponding to different initial conditions as noted before. All other parameters are the same 
as Fig. El 
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Figure 6: The population of all FMO sites at T = 300 K. The initial state is the superposition 
of site 1 and 6. The reorganization energy remains 35 cm -1 . Three different values of phonon 
relaxation time are tested, which are 7" - 1 = 50 fs (a) , 7- 1 = 100 fs (6) and 7" 1 = 166 fs (c). 
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Figure 7: Long time-evolution of the population of each site at T = 300 K, where (a) , (b) and 
(c) correspond to site 1 initially excited, site 6 excited and the superposition of site 1 and 6. The 
reorganization energy is Xj = A = 35 cm -1 , and t" 1 = 7" 1 = 166 fs. 



